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Abstract 

An  Eulerian  and  semi-Lagrangian  finite  element  methods  for  the  solution  of 
the  two  dimensional  advection  equation  were  developed.  Bilinear  rectangular 
elements  were  used.  Linear  stability  analysis  of  the  method  is  given. 

1.  Introduction 

Two  photographs  appearing  in  the  New  York  Times  (March  28,  1994)  show 
the  damage  of  air  pollution  near  big  emission  sources.  But  the  problem  exists 
even  away  from  sources  since  air  pollutants  can  be  transported,  mainly  by  advec¬ 
tion.  Thus  air  pollution  becomes  a  global  problem.  This  physical  phenomenon 
consists  of  three  major  stages  (see  e.g.  Zlatev  [1]): 

1.  emission, 

2.  transport/advection, 

3.  transformation  during  the  transport  which  includes:  diffusion,  deposition 
and  chemical  reactions. 

In  this  paper,  we  only  discuss  the  transport  stage  and  the  solution  of  the 
two  dimensional  advection  equation  by  finite  element  methods. 

2.  Finite  Element  Solution 

The  two  dimensional  advection  equation  is  given  by 

^  ” -^(wc)  -  ^(t;c),  yL<y<yu,  0<t<T  (i) 

at  ax  oy 

where  c  is  the  concentration  of  a  certain  pollutant  and  u  and  v  are  the  wind 
velocity  components  in  the  x  and  y  directions,  respectively.  Clearly,  when  one 
is  interested  in  several  pollutants,  the  equation  is  replaced  by  a  system  of  such 
equations  coupled  only  via  the  chemical  interaction  between  species. 

The  methods  for  numerical  solution  of  the  advection  equation  can  be  divided 
into  five  groups: 

1.  Finite  differences, 

2.  Spectral  methods, 

3.  Finite  volume, 

4.  Characteristic-based  methods  or  semi-Lagrangian, 

5.  Finite  elements. 
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The  finite  difference  methods  are  most  popular  and  have  been  analyzed 
thoroughly  (see  e.g  Richtmeyer  and  Morton  [2]).  Spectral  methods  (see  e.g 
Orszag  [3,4])  are  used  in  weather  forecasting,  but  not  very  much  in  air  pollu¬ 
tion.  The  pseudo  spectral  methods  are  of  the  same  group.  Here  the  solution 
is  approximated  by  a  truncated  polynomial  whose  derivatives  are  substituted 
in  the  equation.  The  spectral  methods  require  periodic  boundary  conditions. 
Finite  volume  or  cell  method  is  based  on  the  integral  form  of  the  equation. 
The  computational  domain  is  divided  into  elements  (volumes  or  cells)  within 
which  the  integration  is  carried  out.  This  method  preserves  the  property  of  con¬ 
servation  (see  Peyret  and  Taylor  [11]).  The  semi-Lagrangian  methods  are  not 
very  popular  among  scientists  working  with  air  pollution  models,  but  these  are 
now  gaining  popularity  in  weather  prediction.  “Discretization  schemes  based 
on  a  semi-Lagrangian  treatment  of  advection  have  elicited  considerable  inter¬ 
est...  since  they  offer  the  promise  of  allowing  larger  time  steps  (with  no  loss  of 
accuracy)  than  Eulerian-based  schemes  whose  time  step  length  is  overly  limited 
by  consideration  of  stability”  (see  Staniforth  and  Cote  [9]).  Semi-Lagrangian 
methods  based  on  finite  difference  or  finite  element  spatial  discretization  were 
developed. 

Here  we  discuss  the  finite  element  approximation  to  the  two  dimensional 
advection  equation.  Both  Eulerian  and  semi-Lagrangian  finite  elements  will  be 
discussed  and  tested.  Software  will  be  available  upon  request  or  electronically 
via  world  wide  web  at  the  URL  address  http://math.nps.navy.mil/--bneta.  The 
advantage  of  finite  elements  is  the  fact  that  the  discretization  can  be  as  easily 
carried  out  for  nonuniform  grids.  Thus  one  can  use  a  fine  grid  only  where  the 
action  is  and  a  coarser  grid  away  from  there.  First  order  linear  one  dimensional 
elements  have  been  previously  used,  see  e.g  Pepper  et  al  [5].  We  now  discuss 
bilinear  finite  elements  on  rectangles.  It  was  shown  by  Neta  and  Williams 
[6]  that  isosceles  triangles  with  linear  basis  functions  and  rectangular  bilinear 
elements  are  superior  to  other  triangulations  and  to  finite  differences.  If  the 
grid  is  uniform,  rectangular  elements  are  preferred  since  Staniforth  et  al  [7] 
have  shown  how  to  evaluate  the  integrals  efficiently  and  the  mciss  matrix  can 
be  replaced  by  a  tensor  product  of  two  tridiagonal  matrices.  If  the  grid  is 
nonuniform  again  the  rectangular  elements  are  preferred,  since  the  isosceles 
triangles  lose  their  shape. 

3.  Bilinear  Finite  Elements 

Discretize  the  rectangular  domain,  by  introducing  the  nodes 
{xi,yj),  i  =  0, 1,...,/+ 1,  i  =  0, 1, . . 1, 

where 

xo  =  XL,  xj+i  =  XR,  2/0  =  VL,  yj+i  =  yu-  (2) 

Suppose  we  number  the  interior  nodes 

(3) 
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from  bottom  left  to  top  right,  see  figure  3  for  the  case  /  =  4,  J  =  3. 
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Figure  3:  node  and  element  numbering 


The  number  of  finite  (rectangular)  elements  is  iVg  =  (/  +  1)(J  +  1)  =  20  in  this 
case. 

We  now  define  the  basis  functions  ipm{^,y)  as  bilinear  functions  on  each 
rectangle,  so  that 


.  .  r  1  at  node  m  . .. 

~  I  Q  other  nodes.  ^  ^ 

To  obtain  the  bilinear  basis  functions  defined  on  the  element,  we  can  make 
a  transformation  of  this  rectangle  to  a  square  centered  at  the  origin  having  sides 
of  length  2  (see  figure  4) . 
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The  transformation  is  given  by 


2  ^ 


_ 2 _  ymfl+yr, 

ym+i-ym^  ym+l~yn 


and  the  basis  functions  in  the  ^  ~  rj  domain  are  given  by 

-!)(??- 1) 

fB  =-5(^  +  l)('?-l) 

(6) 

‘PC  ==  i(^ +  !)(??+ 1) 

PD  = -  1)(»7+ 1) 

where  the  subscripts  denote  the  vertex  at  which  cp  =  1.  Note  that  the  basis 
functions  are  product  of  the  appropriate  linear  basis  functions,  i.e. 

Pa(x,  y)  -  ei(x)em{y) 


= 


?.+i  -  Oi ' 


This  property  is  crucial  to  efficiently  evaluating  the  integrals  (Staniforth  et  al 

[7]). 

The  approximate  problem  becomes 

Mc-Kc=b  (7) 

where  the  entries  of  the  matrices  M,  and  K  are  given  by 


Mij  —  JL  (pjipidxdy 


Kij  — 


The  vector  c  gives  the  concentrations  at  grid  points  at  any  time  t,  and  b  gives 
the  boundary  data 

r  ryu  r^R 

=  /  {U(pi<pj)\l^dy+  {v(pi<pj)\l’[dx  (10) 

i  =  l 

Since  u,  v  are  in  general  functions  of  x  and  y,  we  use  numerical  quadrature  to 
evaluate  Kij,  The  quadrature  we  employed  in  our  case  is  the  two  point  open 
type,  i.e.  ^ 

J  /(2;)da;  =  Y  [/(a  + /i)  + /(a  +  2/i)] ,  (11) 
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where 


(14) 


j_  _  —  xi  u  — 

f^X  —  ^  5  ~  3 

and  S,  the  spacing  for  the  centered  differences  was  arbitrarily  chosen  as 

.Ob{xi+i~xi).  (15) 

Similarly,  we  can  approximate  the  second  integral  in  Kij  except  that  now  the 
points  will  be  (5  =  .05(ym+i  ~  Vm)  units  above  and  below  the  four  points 
E,  F,  G,  H. 

4.  Semi-Lagrangian  Finite  Elements 

Semi-Lagrangian  schemes  belong  to  the  general  class  of  upwinding  methods. 
For  hyperbolic  equations,  upwinding  methods  incorporate  characteristic  infor¬ 
mation  into  the  numerical  method.  In  Lagrangian  schemes,  the  evolution  of  the 
system  is  monitored  by  following  specific  fluid  particles  through  space.  As  a 
result,  Lagrangian  schemes  allow  larger  time  steps  than  Eulerian.  The  problem 
with  fully  Lagrangian  schemes  is  that  an  initially  regularly  spaced  set  of  par¬ 
ticles  will  generally  evolve  into  irregularly  spaced  particles.  As  a  result,  some 
important  features  in  the  flow  may  not  be  captured  properly.  Semi-Lagrangian 
schemes  combine  the  best  of  both  worlds:  the  regular  resolution  of  an  Eulerian 
scheme  and  the  high  stability  of  a  Lagrangian  method.  The  idea  is  to  choose 
a  different  set  of  particles  such  that  at  the  end  of  the  time  step,  they  arrive  at 
points  on  a  regular  Cartesian  grid.  The  departure  points  of  the  particles  are 
determined  by  an  iterative  process  using  the  interpolated  velocity  vector  from 
the  previous  time. 

A  semi-Lagrangian  formulation  of  (1) 

^  [(c«r  +  (cMx  +  CVy)~]  =  0  (16) 

where  is  the  solution  at  the  grid  points  at  time  t  At,  c~  is  the  solution  at 
time  (^- A^)  at  those  points  arriving  at  the  grid  points  at  time  i-(- AL  Since  one 
requires  two  previous  time  levels,  the  program  uses  Matsuno’s  (see  e.g.  Haltiner 
and  Williams  [12])  method  to  get  the  first  time  step. 

In  the  appendix,  we  bring  plots  of  the  solution  for  the  cone  test  (see  e.g. 
Zlatev  [1])  using  Eulerian  finite  elements  with  explicit,  Crank-Nicholson  and 
fully  implicit  time  discretizations  as  well  as  the  semi-Lagrangian  method. 

5.  stability  Analysis 

There  are  four  rectangles  having  a  vertex  in  common,  as  indicated  in  the 
next  figure.  The  approximate  solution  at  the  vertices  of  the  rectangles  may  be 
obtained  by  solving  the  following  first  order  ordinary  differential  equation  (see 
Neta  and  Williams  [6]  for  the  one  dimensional  advection  case). 
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G  H 
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A  B  R 

Figure  6:  rectangular  elements  having  a  common  vertex 


c{P)  +\[6(G)  +  6{B)  +  6{D)  +  i(E)] 

+  1^  [c(-f’)  +  c{H)  +  c(^)  +  c{R)] 

{c(H)  -  c{F)  +  c{R)  -  c{A)  +  4  [c(^)  -  c{D)]} 

-  c(A)  +  4  [c(G)  -  c{B)]}  =  0 

Substitute  a  Fourier  mode 


c(x,y,t)  = 


in  (17)  to  get 


A(<)|  1  +  ^(2  COS  fiAx  +  2  cos  i/Ay)  +  ^4  cos  yAx  cos  z/Ay| 


COS  i/Ay|'2z  sin  fiAx 
+  2  cos  /iAa:|-2?  sin  I'Ay  =  0 


(17) 


(18) 


(19) 
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or 


•  f  u  sin //Ax  V  smi/Ay  \4/.x_ 

(0  +  ^  I  Ax  2  +  cos  //Ax  Ay  2  -i-  cosi/AyJ 


0. 


(20) 


For  the  special  case  of  flow  along  the  x  or  y  axis  one  of  the  terms  in  braces  will 
drop.  For  flow  along  the  diagonal 


u  Ay 
V  Ax 


we  have 


sin  //Ax 


-h  ■ 


sin  i^Ay 


.2  + cos //Ax  2  + cos// Ay 
In  general,  the  ordinary  differential  equation  becomes 

A{t)  +  icrA{t)  =  0, 


}^(<)  = 


0. 


(21) 

(22) 

(23) 


where  cr  is  3  times  the  term  in  braces  in  (20) .  For  the  leap-frog  time  discretiza¬ 
tion 


Aji^i  —  Afi—\  “h  "21(7  At  An  —  0 

(24) 

we  have 

Ai^2  —  —icrAt  ±  >/l  ~  cr2(Af)2, 

(25) 

and  thus  for  stability  (|A|  <  1),  we  must  have 

|crAf|  <  1. 

(26) 

If  we  let 

U  =  max(|i/|,  |t;1) 

(27) 

S  =  min  (Ax,  Ay) 

then  the  method  is  stable  if 

At  2 

T  -W 

(28) 

This  is  the  CFL  condition. 

6.  Fourier  Transform 

The  Fourier  transform  of  a  function  c(x,  y,  t)  is  given  by 

^oo  /*oo 


pOO  pOO 

{c}  c{k,l,t)  =  /  /  c{x,y,t)e~'^'^^'^‘^^dxdy 

«/  — oo  J  —  oo 


(29) 


Taking  the  Fourier  transform  of  the  linearization  of  (1),  one  gets  the  initial  value 
problem 


dc 

dt 


+  i[ku  +  lv)c  —  0, 


(30) 
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(31) 


c{k,L  0)  =  co(A:,/). 

The  solution  of  which  is  given  by 

=  co{kJ)€^‘^\  (32) 


where 

1/ — —(^ku  Iv).  (33) 

To  get  the  solution  c{x,  y,  t)  in  the  physical  domain,  we  have  to  take  the  inverse 
Fourier  transform  and  use  the  convolution  theorem.  This  yields  the  well  known 
solution 

c{x,  y,  t)  =  co{x  -  ut,  y-vt),  (34) 

In  order  to  obtain  the  Fourier  transform  of  the  approximate  solution,  recall 
that 

/oo 

u[x  +  Aa:,  y,  t)e~^^^dx  —  y,  t).  (35) 

-CO 

Applying  Fourier  transform  to  (17),  one  gets 


c  r.  3  sin /: Aa:  .  3  sin /Ay 

- j_  - 1“  IV - 

dt  Aa:  2  +  cos  ^Aar  Ay  2  + cos /Ay 


(36) 


Compare  (36)  and  (30),  to  find  that  k,l  are  replaced  by  (Tx^cTy  respectively, 
where 

3  sin/:Aa:  _  3  sin /Ay 

^  Aa:  2 -f- cos /: Aar  ’  ^  Ay  2-}- cos /Ay 

Note  that  as  Aar  0,  -4  A:  and  as  Ay  -4  0,  (Ty  — )•  /,  thus  at  the  limit  (36) 

becomes  (30).  In  fact 


The  solution  of  (36)  with  the  same  initial  value  is  given  by 

c{k,l,t)  =  CQ{kJ)e^^\  (37) 


where 

i>  =  -[cr^u  +  ayv).  (38) 

The  inverse  transform  is  given  by 

1  rCO  pCO 

c{x,y,t)  =  —  /  (39) 

J  —OO  J—oo 
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or  by  using  convolution 


C[x,y,t)  =  Co  ♦ 


(40) 


7.  Program  Notes 

The  program  can  run  semi-Lagrangian  as  well  as  Eulerian  finite  elements. 
In  the  Eulerian  case  the  time  differencing  is  one  of  the  following: 

•  Explicit 

•  Crank  Nicholson 

•  fully  implicit 

where  the  resulting  linear  system  of  equations  is  solved  by  the  conjugate  gradient 
method  with  symmetric  SOR  preconditioning  (see  e.g.  Ortega  [15]). 

The  four  lines  of  input  contain: 

1.  C  J 

2.  vl,  yu 

3.  AC  T  (final  time  of  integration),  IPLOT  (number  of  time  steps  between 
solution  plots). 

4.  0  (a  parameter  dictating  the  time  integrator).  This  is  needed  only  for 
Eulerian  finite  elements. 

Here  we  include  the  input  file  and  the  programs  used  to  test  the  Eulerian 
and  semi-Lagrangian  finite  element  methods.  First  we  give  the  input  file.  The 
first  line  contains  the  number  of  grid  points  in  the  x  and  y  directions.  The 
second  line  describe  the  rectangular  domain  on  which  the  problem  is  solved. 
The  interval  for  x  is  given  followed  by  the  interval  for  y.  The  third  line  gives 
the  time  step  Af,  the  final  time  of  integration,  and  the  number  of  time  steps 
between  plots.  For  Eulerian  finite  elements,  we  have  a  fourth  line  with  the  value 
of  6^. 

Here  is  the  input  file  for  the  explicit  time  discretization. 

21  21 

0.0  2.0  0.0  2.0 
0.0125  1.0  10 
0. 
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Here  is  the  program  for  the  Eulerian  finite  elements. 

4: - ^ 

*This  program  solves  the  2D  Advection  Equation 

♦  dc/dt  +  d/dx(cu)  +  d/dy(cv)  =  0 
*on  a  square  domain  using  Periodic  B.C's 

+in  both  X  and  y  using  Bilinear  Rectangular  Finite  Elements 
♦and  THETA  Time-Integration  Algorithms. 

♦  THETA=0  ->  FORWARD  EULER  (Explicit) 

♦  THETA=l/2  ->  CRANK-NI COLSON  (Semi-Implicit) 

♦  THETA=2/3  ->  GALERKIN  (Semi-Implicit) 

♦  THETA=1  ->  BACKWARD  EULER  (Implicit) 

♦Written  by  F.X.  Giraldo  on  3/95 

♦  NRC  Fellow 

♦  Department  of  Mathematics 

♦  Naval  Postgraduate  School 

♦  Monterey,  CA  93940 

^ - 1 

program  fem_advect 
implicit  real^8(a-h,o-z) 
parameter  (  imax=21  ) 
c 

c  mxpoi  max  number  of  points 

c  mxele  max  number  of  elements 

c  mxbou  max  number  of  boundary  points 

c  nd  max  number  of  vertices  for  each  elements 

c 

parameter  (  mx=imax*imax,  mxpoi=mx,  mxele=mx,  mxbou=mx/5,  nd=4  ) 
c 

c  global  matrices 
c 

dimension  alhs (mxpoi , mxpoi ) ,  arhs (mxpoi , mxpoi ) ,  b(mxpoi) 
dimension  coord(mxpoi,2) 

integer  intma (mxele ,nd) ,  iboun(mxbou,4) ,  node(imax, imax) 
c 

c  u  velocity  arrays 

c 

dimension  u (mxpoi) 
c 

c  v  velocity  arrays 

c 

dimension  v(mxpoi) 
c 

c  phi  arrays 


14 


c 

dimension  phip(mxpoi),  phiO(mxpoi) 
c 

c  Read  the  Input  Variables  and  create  the  Grid 

c 

c 

c  input  file  contains  4  lines 

c  on  first:  number  of  grid  points  in  x  (nx)  and  y  (ny)  direction 
c  on  second:range  of  x  values  (xmin,  xmax) , 
c  range  of  y  values  (ymin,ymax) 

c  on  third:  delta  t, 

c  final  time  of  integration  (time_f inal) , 

c  number  of  times  steps  between  plots  (iplot) 

c  on  fourth:  theta  (see  above) 
c 

call  init (phiO , u , v , node , coord , intma , iboun , npo in , nelem, 

$  nboun , xmin , xmax , ymin , ymax , ym , nx , ny , nd , dx , dy , dt , 

$  nt  ime ,  theta ,  mxpo  i  ,  mxele ,  mxbou ,  imetx ,  iplot ) 

pi=4. 0*atan(l . 0) 
open ( 1 , f ile=  *  mat lab . out  * ) 
c 

c  isets  =  number  of  time  steps  at  which  solution  is  plotted 
c 

isets=ntime/iplot  +  2 
write( 1 ,*) isets 
c 

c  always  plot  initial  condition 
c 

call  output (phiO,u,v,npoin,t ime, nx,ny,mxpoi) 
c 

c  begin  the  time  marching 

c 

time=0 . 0 
c 

c  CREATE  THE  STIFFNESS  MATRIX  ONCE 

c 

call  Ihs ( alhs , arhs , coord , intma , iboun , node , u , v , npo in , nelem, 

$  nboun , nx , ny , nd , dx , dy , dt , theta , mxpoi , mxele , mxbou , imax ) 

do  i=l,npoin 
do  j=l,npoin 

c  write(2,*)i, j ,alhs(i, j) 

end  do 
end  do 
c 
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c  TIME  MARCH 
c 

do  itime=l ,ntime 
time=time  +  dt 

write(+,^("  timestep  time  =  ",i5,2x,el2.4) Oitime,time/(2.0*pi) 
c 

c  Solve  for  the  GeoPotential 

c 

call  rhs(b,arhs,phiO, coord, intma, iboun,node, 

$  npo in , nelem , nboim , nx , ny , nd , dx , dy , 

$  mxpoi ,mxele ,mxbou, imax) 

if  (theta. eq. 0 . 0)  then 

call  solve_explicit(alhs ,phip,b,npoin, mxpoi) 
else  if  (theta. ne . 0 . 0)  then 

call  pcgm(alhs ,phip,b, npo in, mxpoi) 
end  if 
c 

c  Enforce  B.C.s  Explicitly 

c 

do  i=l,nx 

j l=node(i , 1) 
jny=node(i ,ny) 
phip(jl)=phip(jny) 
end  do 

do  j=l,ny 

il=node(l, j ) 
inx=node(nx, j ) 
ph ip ( i 1 ) =phip ( inx ) 
end  do 
c 

c  Update 

c 

do  ip=l,npoin 

phiO(ip)=phip(ip) 
end  do 
c 

c  check  printing  status 

c 

if  (mod(itime , iplot) . eq. 0) 

$  call  output (phiO ,u,v ,npoin, time ,nx,ny, mxpoi) 

end  do 
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1000  continue 

call  output (phiO , u , v , npoin , t ime , nx , ny ,inxpoi) 

close(l) 

stop 

end 


♦This  subroutine  writes  the  output.  It  is  currently  set  only  to 
♦print  the  concentration  (or  color)  function  at  each  node  point. 
♦Written  by  F.X.  Giraldo  on  2/95 

- - 

subroutine  output (phi ,u , v , npoin , time ,nx ,ny ,mxpoi) 

implicit  real^8(a-h,o-z) 

dimension  phi(mxpoi),  u(mxpoi),  v(mxpoi) 

pi=4. 0^atan(l .0) 

writed,  ^  (2(i6,  lx)  ,el6.8)  O  ,nx,ny,time/(2.0^pi) 
writed,  *  (el2.4)  d  (phi (ip) ,  ip=l, npoin) 
return 
end 


♦This  subroutine  reads  in  the  input  file. 

♦The  info  read  is  the  number  of  grid  points  (in  x  and  y) ,  the  domain, 

♦the  time  step,  the  final  time,  and  the  number  of  time  steps  for  plotting. 
♦Written  by  F.X.  Giraldo  on  2/95 

- - - 

subroutine  init (phiO , uO , vO , node , coord , intma, iboun , npoin , nelem, 

$  nboun,xmin,xmax,ymin,ymax,ym,nx,ny,nd,dx,dy,dt, 

$  nt ime , theta , mxpo i , mxele , mxbou , imax , iplot ) 

implicit  real^8(a-h, o-z) 
dimension  coord (mxpo i, 2) 

dimension  phiO(mxpoi),  uO(mxpoi),  vO(mxpoi) 
integer  intma(mxele ,nd) ,  iboun(mxbou,4) ,  node (imax, imax) 
c 

c  Read  Input  File 

c 

read(*,^)nx,ny 
r ead ( ♦ , ♦ ) xmin , xmax , ymin , ymax 
read ( ♦ , * ) dt , t ime_f inal , iplot 
read(^ ,^)theta 
c 

c  check  bounds 

c 

if  (nx^ny .gt .mxpoi)  then 
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write(*,^(‘'  Error!  -  Need  to  Enlarge  MXP0I")O 
stop 

else  if  ((nx-l)*(ny-l) .gt .mxele)  then 

write(*,^(’'  Error!  -  Need  to  Enlarge  MXELE")0 
stop 

else  if  (2*(nx-l)+2*(ny-l) .gt .mxbou)  then 

write(*,’("  Error!  -  Need  to  Enlarge  MXBOU")’) 
stop 
endif 
c 

c  set  some  constants 
c 

pi=4. 0*ataji(  1 .0) 
nt ime=nint (t ime_f inal/dt ) 
dt=dt*2 . 0*pi 

t ime_f inal=t ime_f inal*2 . 0*pi 
xm=0 . 5* (xmax+xmin) 
ym=0 . 5* (ymax+ymin) 
dx=(xmax-xmin)/ (nx-1) 
dy= (ymax-ymin ) / (ny-1 ) 
xl=xmax-xmin 
yl=ymax-ymin 
w=l  .0 
cx=0 . 25*xl 
cy=0.50*yl 
h=100.0 
rc=0. 125*xl 
velmax=“le5 
c 

c  set  the  Initial  Conditions 
c 

ip=0 

do  j=l,ny 

y=ymin  +  real(j-l)*dy 
do  i=l,nx 

x=xmin  +  real(i-l)*dx 
ip=ip+l 

r=sqrt(  (x“cx)**2  +  (y-cy)**2  ) 
phi0(ip)=0 . 0 
if  (r.lt.rc)  then 

phiO ( ip) =h* ( 1 . 0-r/rc ) 
endif 

u0(ip)=+(y-ym) 

v0(ip)=~(x~xm) 
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vell=u0(ip)**2  +  v0(ip)**2 
velmax=max ( velmax , veil) 
end  do 
end  do 

cfl=dt*sqrt(velmax/(dx**2  +  dy**2)) 
print*, ^  **  CFL  =  *  ,cfl 

GENERATE  COORD 

npoin=nx*ny 

ip=0 

do  j=l,ny 
do  i=l,nx 
ip=ip+l 
node(i, j )=ip 

coordCip, l)=xi  +  dx*real(i-l) 
coord(ip,2)=yi  +  dy*real(j-l) 
end  do 
end  do 

GENERATE  INTMA 

nelem= (nx- 1 ) * (ny- 1 ) 
ie=0 

do  j=l,ny-l 
do  i=l,nx-l 
ie=ie+l 

intma( ie , 1 ) =node ( i , j ) 
intmaC ie , 2 ) =node ( i+1 , j ) 
intma(ie,3)=node(i+l , j+1) 
intma(ie,4)=node(i, j+1) 
end  do 
end  do 

GENERATE  IBOUN 

nboun=2*(nx-l)  +  2*(ny-l) 
ib=0 


bottom  (y=ymin) 
ie=l 

do  i=l,nx-l 
ib=ib+l 

iboun ( ib , 1 ) =node ( i , 1 ) 


iboun(ib,2)=node(i+l , 1) 
iboun(ib,3)=ie 
iboun(ib,4)=l 
ie=ie+l 
end  do 

c  top  (y=ymax) 
ie=nelem 
do  i=nx,2,“l 
ib=ib+l 

iboun(ib, l)=node(i,ny) 
ibonn ( ib , 2 ) =node ( i- 1 , ny ) 
iboun(ib,3)=ie 
iboun(ib,4)=l 
ie=ie‘“l 
end  do 

c  right  (x=xmax) 

ie=(nx-l) 
do  j=l,ny“l 
ib=ib+l 

iboun ( ib , 1 ) =node (nx , j ) 
iboun ( ib , 2 ) =node (nx , j +1 ) 
iboun(ib,3)=ie 
iboun(ib,4)=2 
ie=ie  +  (nx-1) 
end  do 

c  left  (x=xmin) 

ie=nelein  -  (nx-1)  +  1 
do  j=ny,2,“l 
ib=ib+l 

iboun (ib, l)=node(l , j ) 
iboun(ib, 2)=node( 1 , j-l) 
iboun(ib,3)=ie 
iboun (ib, 4) =2 
ie=ie  -  (nx~l) 
end  do 

return 

end 

* - * 

♦This  subroutine  constructs  the  LHS  matrix  for  Bilinear  Rectangular 

♦Elements  for  the  Advection  Equation  with  Periodic 
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♦East-West  and  North-South  Boundary  Conditions. 

♦Written  by  F.X.  Giraldo  on  2/95 

* - 

subrout ine  Ihs (alhs , arhs , coord , intma , iboun , node , u , v , npo in , nelem , 
$  nboun , nx , ny , nd , dx , dy , dt , theta , mxpo i , mxele , mxbou , imax ) 

implicit  real^8(a-h,o-z) 
parameter  (mx=3000) 
c 

c  global  arrays 
c 

dimension  alhs(mxpoi ,mxpoi) ,  arhs (mxpo i , mxpo i) 
dimension  coord(mxpoi ,2) ,  u(mxpoi),  v(mxpoi) 
integer  intma(mxele ,nd) ,  node (imax, imax) ,  iboun (mxbou, 4) 
c 

c  local  coordinate  system  for  a  CCW  ordered  rectangle 
c 

dimension  xi(4),  eta(4) ,  a_temp(mx) 
data  xi  /“I,  1,  1,”1  / 
data  eta  1,1/ 

c 

c  initialize  the  global  matrix 
c 

do  j=l,npoin 
do  i=l,npoin 

alhs(i, j )=0.0 
arhs ( i, j ) =0.0 
end  do 
end  do 


c 

c 

c 

c 

c 

c 


loop  thru  the  elements 


do  ie=l, nelem 


assemble  element  matrix  cm  ==  consistent  mass  matrix 

ckc  ==  conduction-like  matrix 


do  i=l,nd 

ii=intma(ie, i) 
cm_lump=0 . 0 
do  j=l,nd 

j j=intma(ie , j ) 

cm=(2,0+2.0/3.0^xi(i)^xi(j) )^(2.0+2.0/3.0^eta(i)^eta(j )) 
ckx=0 . 0 
cky=0 . 0 
do  k=l,nd 


21 


uk=u(intma(ie ,k) ) 
vk=v(intma(ie ,k) ) 
ckx=ckx  +  uk* 

$  (2.0*xi(i)  +  2.0/3.0*xi(i)*xi(j)+xi(k))* 

$  (2 . 0+2 . 0/3 . 0*(eta(i)*eta( j )  +  eta(i)*eta(k)  +  eta( j ) ♦eta(k) ) ) 

cky=cky  +  vk* 

$  (2.0*eta(i)  +  2 . 0/3 . 0*eta(i)*eta( j ) ♦eta(k) )* 

$  (2.0+2.0/3.0*(xi(i)*xi(j)  +  xi(i)*xi(k)  +  xi( j ) ♦xi(k) ) ) 

end  do 

cinin=dx*dy/ 64 . 0*cm 

ck=dy/128 . 0*ckx  +  dx/128.0*cky 

alhs(ii , j j )=alhs(ii, j j )  +  l.O/dt+cmm  -  theta*ck 
arhsCii,  j  j  )=arhs(ii,  j  j  )  +  1.0/dt*cmm  +  (1 . 0’~theta)+ck 
end  do 
end  do 
end  do 
c 

c  Account  for  Periodicity 
c 

do  i=l,nx 

j l=node(i , 1) 
jny=node(i,ny) 
do  jj=l,npoin 

a_teiBp(jj)=alhs(jl,jj)  +  alhs(jny,jj) 
end  do 

do  jj=l,npoin 

alhsCjl, jj)=a_temp(jj) 
alhs(jny, jj)=a_temp(jj) 
end  do 

alhsCjl, jl)=a_temp(jl)  +  a_temp(jny) 
alhsCjny , jny)=a„temp(j 1)  +  a_temp(jny) 
end  do 

do  j=l,ny 

il=node(l,j) 
inx=fiode(nx ,  j ) 
do  jj=l,npoin 

a_temp(j j )=alhs(il, j j )  +  alhs(inx,jj) 
end  do 

do  jj=l,npoin 

alhsCil, jj)=a_temp(jj) 
alhs(inx, j j )=a_temp( j j ) 
end  do 

alhs (il , il)=a_temp(il)  +  a_temp(inx) 
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alhs(inx,inx)=a^temp(il)  +  a^temp(inx) 
end  do 
c 

c  If  theta=0,  then  Lump, 
c 

if  (theta. eq. 0 . 0)  then 
do  i=l,npoin 
asum=0 . 0 
do  j=l,npoin 

asum=asum  +  alhs(i,j) 
alhs(i, j)=0.0 
end  do 

alhs(i,i)=asum 
end  do 
endif 

return 

end 


♦The  next  4  Subroutines  solve  a  linear  system  [A]{x}-{b}',  for  n  unknowns  using 
♦the  CONJUGATE  GRADIENT  METHOD  with  an  SSOR  preconditioner. 

♦The  variable  w  determines  the  preconditioner,  for  w=l  it  is 
♦Symmetric  Jacobi  and  for  w>l  it  is  Symmetric  SOR  (SSOR). 

♦Written  by  F.X.  Giraldo  on  4/14/90 
♦Modified  for  any  type  of  spd  Matrix  A  on  2/95 


subroutine  pcgm(a,x,b,n,mx) 

implicit  real^8(a-h,o-z) 

parameter  (nmax=3000,  kmax=15,  w=1.0) 

dimension  a(mx,mx),  b(mx),  x(mx) 

dimension  r(nmax) ,  rw(nmax) ,  p(nmax),  ap(nmax) 

do  i=l,n 
sum=0 . 0 
do  j=l,n 

STim=sum  +  a(i,j)^x(j) 
end  do 

r(i)=b(i)  -  sum 
end  do 

call  ssor(r ,rw,a,n,w,nmax,mx) 

do  i=l,n 

p(i)=rw(i) 
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end  do 


call  ip(rop,rw,r ,n,nmax) 

do  k=l,kmax 

do  i=l,n 
sum=0 . 0 
do  j=l,n 

sum=sum  +  a(i,j)*p(j) 
end  do 
ap(i)=sum 
end  do 

alfden=0 . 0 
do  i=l,n 

alfden=alfden  +  p(i)*ap(i) 
end  do 

alf=-rop/alfden 
do  i=l,n 

x(i)=x(i) -alf  *p ( i ) 
end  do 

do  i=l,n 

r(i)=r(i)  +  alf*ap(i) 
end  do 

call  convtest (rop,r ,n,f lag,ninax) 
if  (flag  .eq.  1.0)  goto  200 

call  ssor(r ,rw,a,n,w,nmax,mx) 
call  ip(rnp,rw,r ,n,nmax) 

beta=rnp/rop 
do  i=l,n 

p(i)=rw(i)  +  beta*p(i) 
end  do 
rop=rnp 
end  do 
k=k-l 

print*,*  No  Convergence* 

200  continue 
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print*,’  SSOR  PCG  Loops  =  ’ ,k 

return 

end 


subroutine  ssor(r , rw, a,n, w,ninax,mx) 
implicit  real*8(a-h, o-z) 
dimension  r(nmax) ,  rw(nmax) 
dimension  a(mx,mx) 

*  symmetric  sor  on  the  residual 

*  zeroing  the  wiggle  residual 

do  i=l,n 

rw(i)=0,0 
end  do 

*  up  sweep 

do  i=l,n 
sum=0 . 0 
do  j=l,n 

sum=sum  +  a(i,j)*rw(j) 
end  do 

rw(i)=rw(i)  +  w/a(i,i)*(  r(i)  -  sum  ) 
end  do 

*  down  sweep 

do  i=n,l,~l 
sum=0 . 0 
do  j=l,n 

sum=sum  +  a(i,j)*rw(j) 
end  do 

rw(i)=rw(i)  +  w/a(i,i)*(  r(i)  -  sum  ) 
end  do 
return 
end 


subroutine  ip(f ,fw,f r ,n,nmax) 
implicit  real+SCa-'hjO-z) 
dimension  fw(nmajc) ,  fr(nmax) 
f=0.0 
do  i=l,n 

f=f  +  fw(i)*fr(i) 
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end  do 
return 
end 

- ♦ 

subrout ine  convtest (rop , r , n , f lag ,nmax) 
implicit  real*8(a-h,o-z) 
dimension  r(nmax) 

emin=l . Oe-6 
f  lag==0 . 0 
rtest=0 . 0 

if  (rop  .It.  emin)  then 
do  i=l,n 

rtest=rtest  +  r(i)**2.0 
end  do 

if  (rtest  .It.  emin)  flag=1.0 
endif 
return 
end 

* - * 

♦This  subroutine  Builds  the  RHS  vector  for  Bilinear  Rectangular  Finite 
♦Elements  for  the  2D  SLSI  Shallow  Water  Equations  with  Periodic 
♦West-East  and  North-South  Boundaries . 


♦Written  by  F.X.  Giraldo  on  2/95 

^ - * 

subrout ine  rhs (b , arhs , phiO , coord , intma , iboun , node , 

$  npo in , nelem , nboun , nx , ny , nd , dx , dy , 

$  mxpoi,mxele,mxbou,imax) 


implicit  realms (a-h,o-z) 
parameter  (mx=3000) 
c 

c  global  arrays 

c 

dimension  b(mxpoi),  arhs (mxpoi ,mxpoi) ,  phiO(mxpoi) 
dimension  coord(mxpoi ,2) 

integer- intma(mxele,nd) ,  iboun(mxbou,4) ,  node(imax, imax) 
c 

c  local  coordinates  for  a  CCW  ordered  Rectangle 

c 

dimension  xi (4) ,  eta(4) 
data  xi  /-I,  1,  1,-1/ 
data  eta  /-1,-1,  1,  1/ 
c 

c  initialize  the  right  hand  side 
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c 


do  ip=l,npoin 
b(ip)=0 . 0 

end  do 
c 

c  Construct  the  RHS  Vector 

c 

do  i=l,npoin 
do  j=l,npoin 

b(i)=b(i)  +  arhs(i, j)*phiO(j) 
end  do 
end  do 
c 

c  Account  for  Periodicity 

c 

do  i=l,nx 

j l=node(i , 1) 
jny=node(i ,ny) 
b_temp=b( j 1)  +  b(jny) 
b(j l)=b_temp 
b( jny)=b_temp 
end  do 
do  j=l,ny 

il=node(l , j ) 
inx=node(nx, j ) 
b_temp=b(il)  +  b(inx) 
b(il)=b_temp 
b(inx)=b_tenip 
end  do 
return 
end 

- - - 

♦This  subroutine  solves  a  Linear  NxN  system  where  the  Coefficient 
♦Matrix  is  diagonal. 

♦Written  by  F.X.  Giraldo  on  4/14/90 

* - - - - - ♦ 

subroutine  solve_explicit (a,x,b,n,mx) 
implicit  realms (a-h,o-z) 
dimension  a(mx,mx),  b(mx) ,  x(mx) 
do  i=l,n 

x(i)=b(i)/a(i,i) 
end  do 
return 
end 
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Here  is  the  program  for  the  semi-Lagrangian  finite  elements. 

♦ - * 

♦This  program  solves  the  2D  Advection  Equation 
♦  dc/dt  +  d/dx(cu)  +  d/dy(cv)  =  0 

♦on  a  square  domain  using  Periodic  B.C.'s 
♦in  both  X  and  y  and  using  Semi-Implicit  Semi-Lagrangian 
♦Bilinear  Rectangular  Finite  Elements. 

♦Written  by  F.X.  Giraldo  on  3/95 


♦  NRC  Fellow 

♦  Department  of  Mathematics 

♦  Naval  Postgraduate  School 

♦  Monterey,  CA  93940 

♦  - ♦ 


program  slt_advect 
implicit  realms (a-h, o-z) 
parameter  (  imax=21  ) 
c 

c  mxpoi  max  number  of  points 

c  mxele  max  number  of  elements 

c  mxbou  max  number  of  boundary  points 

c  nd  max  number  of  vertices  for  each  elements 

c 

parameter  (  mx=imax^imax,  mxpoi=mx,  mxele=mx,  mxbou=mx/5,  nd=4  ) 
c 

c  global  matrices 
c 

dimension  a(mxpoi , mxpoi) ,  b(mxpoi) 
dimension  f (mxpoi) 
dimension  coord(mxpoi , 2) 
dimension  cmat (mxpoi) 

integer  intma(mxele ,nd) ,  iboun(mxbou,4) 
c 

c  u  velocity  arrays 

c 

dimension  urn (mxpoi ) ,  uO(mxpoi),  up(mxpoi) 

dimension  u0_x (mxpoi) 

c 

c  V  velocity  arrays 

c 

dimension  vm(mxpoi) ,  vO(mxpoi),  vp(mxpoi) 

dimension  v0_y (mxpoi) 

c 

c  phi  arrays 
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c 

dimension  phim(mxpoi) ,  phiO(mxpoi),  phip(mxpoi) 
c 

c  departure  point  arrays 
c 

dimension  alfm(mxpoi ,2) ,  alf0(mxpoi,2) 
c 

c  spline  derivative  arrays 

c 

dimension  dphim(mxpoi) 
dimension  duO(mxpoi),  duO_x(mxpoi) 
dimension  dvO(mxpoi) ,  dvO_y(mxpoi) 
c 

c  Auxiliary  Matrices 

c 

int  eger  node ( imax , imax ) 
c 

c  Read  the  Input  Variables  and  create  the  Grid 

c 

c  input  file  contains  4  lines 

c  on  first:  number  of  grid  points  in  x  (nx)  and  y  (ny)  direction 
c  on  second:range  of  x  values  (xmin,  xmax), 
c  range  of  y  values  (ymin,ymax) 

c  on  third:  delta  t, 

c  final  time  of  integration  (time_f inal) , 

c  number  of  times  steps  between  plots  (iplot) 

c 

call  init (phiO , uO , vO , node , coord , alf 0 , intma , iboun , npo in , 

$  nelem , nboun , xmin , xmax , ymin , ymax , nx , ny , dx , dy , dt , 

$  nt ime , ym , nd , mxpo i , mxele , mxbou , imax , iplot ) 

c 

c  Construct  the  Lumped  Mass  Matrix  used  for  the  derivative  computation 
c 

call  get _geom(nelem, npo in, intma, coord, cmat ,node, 

$  nx,ny,nd,dx,dy,mxpoi, mxele, imax) 

time=0 . 0 
pi=4.0*atan(1.0) 
open ( 1 , f ile=  *  mat lab . out  * ) 
c 

c  isets  =  number  of  time  steps  at  which  solution  is  plotted 
c 

isets=ntime/iplot  +  2 
write(l,*)isets 


c 
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always  plot  initial  condition 

call  outputCphiO ,u0 , vO ,npoin,tinie,nx,ny ,mxpoi) 
begin  the  time  marching 

Do  the  1st  Time-step  Integration  to  get  2-time  levels 

do  itime=l,l 

time=time  +  dt 
dtime=time/(2 . 0*pi) 

write(*,^("  timestep  time  =  " , i5,2x,el2.4) O itime,dtime 
call  matsuno(phim,phiO ,phip,um,uO ,iip, 

$  vm, vO,vp, coord, intma,npo in, nelem,dt ,dx,dy, 

$  node,nx,ny ,ym,nd,mxpoi,mxele, imax) 

call  update (phim , phiO , phip , urn , uO , up , vm , vO , vp , alf m , alf 0 , 
$  npoin,mxpoi) 

end  do 


TIME  MARCH 

do  itime=2 ,ntime 

time=time  +  dt 
dtime=time/(2 . 0*pi) 

write(*,’("  timestep  time  =  " , i5,2x,el2.4) Oitime,dtime 
1st,  DETERMINE  DEPARTURE  POINT 

call  splie2(du0,u0,coord,nx,ny,node,mxpoi, imax) 
call  splie2(dv0, vO, coord, nx,ny, node, mxpoi, imax) 
call  depart (alf 0 , alf m, coord , uO , duO , vO , dvO , 
i  npoin,dt ,xmin,xmax,ymin,yraax, 

i  nd,mxpoi,node,nx,ny , imax) 

2nd , *  COMPUTE  DERIVATIVES 

call  der iv_x (uO_x , uO , intma , cmat , node , npoin , nelem , 
i  nx,ny,nd,dy , mxpoi, mxele, imax) 

call  deriv_y(vO_y , vO , intma, cmat , node, npoin, nelem, 
i  nx,ny,nd,dx, mxpoi, mxele, imax) 


3rd,  INTERPOLATE  PHIM,  UM^X=U0_X,  VM_Y=V0^Y 
at  the  departure  point  =  X  -  2* ALPHA 
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♦ 

call  splie2 (dphim , phim , coord , nx , ny ,node ,mxpoi , imax) 
call  splie2(du0_x , uO_x , coord , nx , ny ,node ,mxpoi , imax) 
call  splie2 (dvO_y , vO^y , coord , nx , ny ,node ,mxpoi , imax) 
c 

c  INTERPOLATE  THE  RIGHT  HAND  SIDE  FUNCTION 

c 

call  int erp (f , phim , uO_x , vO_y , dphim , duO_x , dvO_y , coord , 

$  alf  0 , node , npoin , dt , xmin , xmax , ymin , ymax , yra , 

$  nx,ny ,nd,mxpoi , imax) 

do  ip=l, npoin 
phip(ip)=f (ip) 
up(ip)=uO(ip) 
vp(ip)=vO(ip) 
end  do 
c 

c  Update  the  values 

c 

call  update (phim , phiO , phip , urn ,u0 , up , vm , vO , vp , alf m , alf 0 , 

$  npoin,mxpoi) 

c 

c  check  time  for  printing  output 
c 

if  (mod(itime , iplot) .eq. 0) 

$  call  output ( ph iO , uO , vO , npo in , t ime , nx , ny , mxpo i ) 

end  do 

1000  continue 

call  output (phiO , uO , vO , npo in , t ime , nx , ny , mxpo i ) 

close(l) 

stop 

end 

- - 

♦This  subroutine  finds  the  departure  point 

+  ALPHA 1=DT*U ( X-ALPHAl , Y-ALPHA2 , T )  ALPHA2=DT*V (X-ALPHAl , Y- ALPHA2 , T ) 
♦Written  by  F.X.  Giraldo  on  2/95 

- - 

subrout ine  depart ( alf 0 , alf m , coord , uO , duO , vO , dvO , 

$  npo in, dt , xmin, xmax, ymin, ymax, 

$  nd , mxpo i , node , nx , ny , imax ) 

implicit  real*8(a-h,o-'z) 
dimension  alf 0(mxpoi , 2) ,  alfm(mxpoi ,2) 


♦ 


* 


31 


dimension  uO(mxpoi),  duO(mxpoi) 
dimension  vO(mxpoi),  dvO(mxpoi) 
dimension  coord(mxpoi , 2) 
integer  node(imax , imax) 


do  ip=l,npoin 

alphal=alfm(ip , 1) 
alpha2=alfra(ip ,2) 
do  k=l,3 

xd=coord(ip, 1)  -  alphal 
yd=coord(ip,2)  -  alpha2 
if  (xd.lt.xmin)  xd=xmax  “ 
if  (xd.gt.xmax)  xd=xmin  + 
if  (yd.lt.ymin)  yd=yraax  - 
if  (yd.gt.ymax)  yd=ymin  + 


(xmin  -  xd) 
(xd  -  xmax) 
(ymin  -  yd) 
(yd  -  ymax) 


if  (  (xd. It .xmin. or .xd.gt .xmax)  .or. 

$  (yd.lt .ymin. or. yd. gt .ymax)  )  then 

print*,*  Error  in  DEPART* 
print*, *XD  out  of  Range  =  *,xd,yd 
print*,* ip  alpha  =  *, ip, alphal , alpha2 
stop 
endif 


cal 1  spl in2 ( u , coord , uO , duO , xd , yd , nx , ny , node , mxpo i , imax ) 
call  spl in2 ( V, coord, vO,dvO ,xd, yd, nx,ny, node , mxpo i , imax) 

alphal=dt*u 
alpha2=dt*v 
end  do 

alfO(ip, l)=alphal 
alf 0 ( ip , 2 ) =alpha2 
end  do 


return 

end 

+ - 1 

*This  subroutine  computes  the  4th  Order  Accurate  derivative  WRT  X  of  the 
*variable  UNKNO  and  stores  it  in  DERIP  using  Bilinear  Rectangular 
*Finite  Elements 

*Written  by  F.X.  Giraldo  on  2/95 

- ♦ 

subroutine  deriv_x (derip , unkno , intma , cmat , node , npoin ,nelem , 

$  nx,ny ,nd,dy, mxpo i ,mxele, imax) 
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implicit  real*8(a-h,o-z) 
c 

c  global  arrays 
c 

dimension  derip (mxpoi) ,  unkno(mxpoi) 
dimension  cmat (mxpoi) 

integer  intma(mxele,nd) ,  node ( imax , imax ) 
c 

c  local  coordinate  system  for  a  CCW  ordered  rectangle 
c 

dimension  xi(4),  eta(4) 
data  xi  /  **1,  1,  1,-1  / 
data  eta  /  -1,-1,  1,1/ 
c 

c  initialize  arrays 

c 

do  ip=l,npoin 
derip(ip)=0 . 0 
end  do 
c 

c  Loop  thru  the  Elements  and  find  the  ISt  DERIVATIVES 
c 

do  ie=l,nelem 
do  i=l,nd 
phi_x=0 . 0 
ip=intma(ie, i) 
do  k=l,nd 

phik=unkno(intma(ie,k) ) 

phi^x=phi_x  +  xi(k)*phik*(  2.0  +  2.0/3.0*eta(k)*eta(i)  ) 
end  do 

derip(ip)=derip(ip)  +  dy/16 . 0*phi_x 
end  do 
end  do 
c 

c  Account  for  Periodicity 
c 

do  j=l,ny 

il=node(l, j) 
inx=node(nx, j ) 

derip_temp=derip(il)  +  derip(inx) 
derip ( i 1 ) =der ip_temp 
derip (inx)=derip_temp 
end  do 
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c  Now  multiply  by  the  Inverse  of  the  Lumped  mass  matrix 
c 

do  ip=l,npoin 

derip(ip)=derip(ip)*cmat (ip) 
end  do 

return 

end 

* - * 

★This  subroutine  computes  the  4th  Order  Accurate  derivative  WRT  Y  of  the 
♦variable  UNKNO  and  stores  it  in  DERIP  using  Bilinear  Rectangular 
♦Finite  Elements 

♦Written  by  F.X.  Giraldo  on  2/95 

^ - ^ 

subrout ine  deriv_y (derip , unkno , intma , cmat , node , npo in , nelem , 

$  nx,ny,nd,dx,mxpoi,mxele,imax) 

implicit  real^8(a-h,o-z) 
c 

c  global  arrays 
c 

dimension  derip(mxpoi) ,  unkno(mxpoi) 
dimension  cmat(mxpoi) 

integer  intma(mxele ,nd) ,  node(imax, imax) 
c 

c  local  coordinate  system  for  a  CCW  ordered  rectangle 
c 

dimension  xi(4),  eta(4) 
data  xi  /  -1,  1,  1,~1  / 
data  eta  /  -1,-1,  1,1/ 
c 

c  initialize  arrays 
c 

do  ip=l,npoin 
derip(ip)=0 . 0 
end  do 
c 

c  Loop  thru  the  Elements  and  find  the  ISt  DERIVATIVES 
c 

do  ie=l, nelem 
do  i=l,nd 
phi_y=0 . 0 
ip=intma(ie , i) 
do  k=l,nd 

phik=unkno(intma(ie,k) ) 
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phi_y=phi_y  +  eta(k) ♦phik* (  2.0  +  2 . 0/3 . 0*xi(k)*xi(i)  ) 
end  do 

derip(ip)=derip(ip)  +  dx/16 . 0*phi_y 
end  do 
end  do 
c 

c  Account  for  Periodicity 
c 

do  j=l,ny 

il=node(l, j ) 
inx=node(nx, j ) 

derip_temp=derip(il)  +  derip(inx) 
derip ( i 1 ) =der ip_t emp 
derip (inx)=derip_temp 
end  do 
c 

c  Now  multiply  by  the  Inverse  of  the  Lumped  mass  matrix 
c 

do  ip=l,npoin 

derip(ip)=derip(ip)*cmat(ip) 
end  do 

return 

end 


♦This  subroutine  computes  the  Inverse  Lumped  Mass  Matrix 

♦for  Bilinear  Rectangular  Finite  Elements  used  for  obtaining  the 

♦4th  Order  Accurate  1st  and  2nd  derivatives  in  both  X  and  Y. 

♦Written  by  F.X.  Giraldo  on  2/95 

♦ - - 

subrout ine  get_geom(nelem , npoin , intma , coord , cmat , node , 

$  nx , ny , nd , dx , dy , mxpo i , mxe le , imax ) 

implicit  real*8(a-h,o-z) 
c 

c  global  arrays 

c 

dimension  coord(mxpoi ,2) ,  cmat(mxpoi) 
integer  intma(mxele ,nd) ,  node (imax, imax) 
c 

c  Compute  the  inverse  lumped  mass  matrix 

c 

do  ip=l, npoin 
cmat (ip) =0.0 
end  do 
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do  ie=l,nelem 
do  in=l,nd 

ip=intma(ie , in) 

Croat (ip)=cmat (ip)  +  dx*dy/4.0 
end  do 
end  do 
c 

c  Account  for  Periodicity  in  X 
c 

do  j=l,ny 

il=node(l , j ) 
inx=node(nx, j ) 

cmat_temp=cmat (il)  +  Croat (inx) 

Croat ( i 1 ) =cmat  _t  erop 
Croat ( inx ) =croat_temp 
end  do 
c 

c  Invert  the  lumped  mass  matrix 
c 

do  ip=l,npoin 

cmat (ip)=l . 0/cmat (ip) 
end  do 

return 

end 

* - * 

♦This  subroutine  reads  in  the  input  file. 

♦The  info  read  is  the  number  of  grid  points  (in  x  and  y) ,  the  domain, 

♦the  time  step,  the  final  time,  coid  the  number  of  time  steps  for  plotting. 


♦Written  by  F.X.  Giraldo  on  2/95 

* - ^ 

subroutine  init (phiO , uO , vO ,node , coord, alfO , intma, iboun ,npoin , 

$  nelem , nboun , xmin , xmax , ymin , ymax , nx , ny , dx , dy , dt , 

$  ntime ,ym,nd,mxpoi,mxele ,mxbou, imax, iplot) 


implicit  real^8(a-h,o-z) 
dimension  coord(mxpoi ,2) ,  alf 0(mxpoi,2) 
dimension  phiO(mxpoi),  uO(mxpoi),  vO(mxpoi) 
integer  intma(mxele,nd) ,  iboun(mxbou,4) ,  node (imax, imax) 
c 

c  Read  Input  File 
c 

read(^ , ♦)nx,ny 

read ( ♦ , ♦ ) xmin , xmax , ymin , ymax 
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read ( * , * ) dt , t ime^f inal , iplot 
c 

c  check  bounds 
c 

if  (nx*ny ,gt .mxpoi)  then 

write(*,^("  Error!  -  Need  to  Enlarge  MXP0I’')O 
stop 

else  if  ((nx-l)*(ny-“l)  .gt  .mxele)  then 

write(*,^("  Error!  -  Need  to  Enlarge  MXELE") O 
stop 

else  if  (2*(nx-l)+2*(ny--l)  .gt  .mxbou)  then 

write(*,'("  Error!  -  Need  to  Enlarge  MXB0U")O 
stop 
endif 
c 

c  set  some  constants 
c 

pi=4.0*atan(l .0) 

nt  ime=nint  ( t  ime__f  inal/ dt ) 

dt=dt*2.0*pi 

time_f inal=time_f inal*2 .0*pi 

xm=0 . 5* (xmax+xmin) 

ym=0 . 5* (ymax+ymin) 

dx=(xmax~xmin)/ (nx-l) 

dy=(ymax”ymin)/ (ny-1) 

xl=xmax“-xmin 

yl=ymax-ymin 

w=1.0 

cx=0.25*xl 
cy=0.50*yl 
h=100.0 
rc=0 . 125*xl 
velmax=-le5 
c 

c  set  the  Initial  Conditions 
c 

ip=0 

do  j=l,ny 

y=ymin  +  real(j-l)*dy 
do  i=l,nx 

x=xmin  +  real(i-l)*dx 
ip=ip+l 

r=sqrt(  (x”cx)**2  +  (y-cy)**2  ) 
phi0(ip)=0 . 0 
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if  (r.lt.rc)  then 

phiO(ip)=h*(l . 0-r/rc) 
endif 

uO(ip)=+(y-yin) 
vO(ip)=-(x-xm) 
vell=u0(ip)**2  +  v0(ip)**2 
velmax=inax  (  velmax ,  vel  1 ) 
end  do 
end  do 

cf l=dt*sqrt (velraax/(dx**2  +  dy**2)) 
print*,  ^  **  CFL  =  \cfl 

GENERATE  COORD 

npoin=nx*ny 

ip=0 

do  j=l,ny 
do  i=l,nx 
ip=ip+l 
node(i, j )=ip 

coord(ip, l)=xi  +  dx*real(i~l) 
coord(ip,2)=yi  +  dy*real(j-l) 
alf0(ip,l)=0.5*dx 
alf  0 ( ip , 2 ) =0 . 5*dy 
end  do 
end  do 

GENERATE  INTMA 

nelem=(nx“l)*(ny-l) 

ie=0 

do  j  =  l,ny-“l 
do  i=l,nx“l 
ie=ie+l 

intma( ie , 1 ) =node ( i , j ) 
intma ( i e , 2 ) =node ( i+ 1 , j ) 
intma(ie,3)=node(i+l , j+1) 
intma(ie,4)=node(i , j+1) 
end  do 
end  do 

GENERATE  IBOUN 


nboun=2* (nx^ 1 ) 


ib=0 


c  bottom  (y=ymin) 

ie=l 

do  i=l,nx”l 
ib=ib+l 

iboun(ib, l)=node(i , 1) 
iboun(ib,2)=node(i+l ,1) 
iboun(ib,3)=ie 
iboun(ib,4)=3 
ie=ie+l 
end  do 

c  top  (y=ymax) 

ielem=nelem 
do  i=nx,2,-‘l 
ib=ib+l 

iboun(ib, l)=node(i ,ny) 
iboun(ib,2)=node(i-l ,ny) 
iboun(ib,3)=ielem 
iboun(ib,4)=3 
ielem=ielem-l 
end  do 

return 

end 


♦This  subroutine  interpolates  the  right  hand  side  function  for  the 
♦2D  Advection  Equation  using  a  3-time  level  Semi-Lagrangian 
♦Semi-Implicit  Method 
♦Written  by  F,X.  Giraldo  on  2/95 

♦ - 

subrout ine  int erp (f , phim , uO_x , vO_y , dphim , duO^x , dvO^y , coord , 

$  alf 0 , node , npo in , dt , xmin , xmax , ymin , ymax ,  ym , 

$  nx , ny , nd , mxpo i , imax ) 

implicit  realms (a-h,o-z) 
dimension  f(mxpoi) 

dimension  phim(mxpoi),  uO_x(mxpoi),  vO_y(mxpoi) 
dimension  dphim(mxpoi) ,duO_x(mxpoi) ,dvO_y(mxpoi) 
dimension  coord(mxpoi,2) ,  alf 0(mxpoi ,2) 
integer  node (imax, imax) 

do  ip=l,npoin 
c 
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c  1st,  Interpolate  values  =  F(  x  -  2*alpha,  t  -  dt  ) 
c 

xd=coordCip, 1)  -  2 . 0*alf0(ip, 1) 
yd=coord(ip,2)  -  2 . 0*alf 0(ip, 2) 
if  (xd . It . xmin)  xd=xmax  -  (xmin  -  xd) 
if  (xd.gt  .xmax)  xd=xinin  +  (xd  -  xmax) 
if  (yd.lt.ymin)  yd=ymax  -  (ymin  -  yd) 
if  (yd.gt . ymax)  yd=ymin  +  (yd  -  ymax) 

if  (  (xd. It .xmin. or .xd.gt .xmax)  .or. 

$  (yd.lt.ymin.or .yd.gt .ymax)  )  then 

print*, ^  Error  in  INTERP' 
print*, ^XD  out  of  Range  =  ',xd,yd 
stop 
endif 

c  Do  Interpolation 

call  spl in2 (phi, coord, phim,dphim,xd, yd, nx,ny, node, mxpoi, imaix) 
call  splin2(u_x, coord, uO_x,duO_x,xd, yd, nx,ny, node , mxpoi, imax) 
call  splin2 ( v_y , coord , vO_y , dvO_y , xd , yd , nx , ny , node , mxpoi , imajc) 

f (ip)=( 1 . 0-dt* (u_x  +  v_y) )/( 1 . 0+dt*(u0_x(ip)  +  vO_y(ip) ) )*phi 
end  do 

return 

end 

* - :fc 

*This  subroutine  solves  the  2D  Advection  Equation 

*using  the  Backward  Euler  with  a  predictor-corrector  strategy 


*Written  by  F.X.  Giraldo  on  2/95 

^ - ^ 

subrout ine  mat suno (phim , phiO , phip , urn , uO , up , 

$  vm, vO ,vp, coord, intma,npoin,nelem,dt ,dx,dy, 

$  node, nx,ny,ym,nd, mxpoi ,mxele, imax) 


implicit  real*8(a-h, o-z) 

dimension  ph im (mxpoi) ,  phiO(mxpoi),  phip(mxpoi) 
dimension  um (mxpoi ) ,  uO(mxpoi),  up(mxpoi) 
dimension  vm(mxpoi),  vO(mxpoi),  vp(mxpoi) 
dimension  coord(mxpoi , 2) 
integer  intma(mxele ,nd) ,  node (imax, imax) 
c 

c  Loop  through  the  points  and  integrate  using  Forward  Time 
c  and  Centered  Space. . . 
c 
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Predictor  Stage  (forward  Euler) 
do  j=l,ny 

j2=j+l 

if  (jl.lt. 1)  jl=ny-l 
if  (j2.gt.ny)  j2=2 
do  i=l,nx 
il=i-l 
i2=i+l 

if  (il.lt. 1)  il=nx-l 
if  (i2.gt.nx)  i2=2 

Set  up  the  nodes  in  X  and  Y 

ip=node(i , j ) 
ipl=node(il, j ) 
ip2=node(i2, j ) 
jpl=node(i , j 1) 
jp2=node(i, j2) 

integrate  PHI 

ph  im  (  ip )  =ph  iO  (  ip  ) 

$  -0.5*dt/dx*u0(ip)*(  phi0(ip2)-phi0(ipl)  ) 
$  -0.5*dt/dy*v0(ip)*(  phiO( jp2)-phi0(jpl)  ) 

integrate  U 

um(ip)=uO(ip) 

integrate  V 

vm(ip)=vO(ip) 
end  do 

end  do 


Corrector  Stage  (backward  Euler) 


if  (jl.lt. 1)  jl=ny-l 
if  (j2.gt.ny)  j2=2 
do  i=l,nx 
il=i-l 
i2=i+l 

if  (il.lt.  1)  il=nx-l 
if  (i2.gt.nx)  i2=2 
c 

c  Set  up  the  nodes  in  X  and  Y 

c 

ip=node(i, j) 
ipl=node(il , j ) 
ip2=node(i2, j ) 
jpl=node(i, j 1) 
jp2=node(i , j2) 
c 

c  integrate  PHI 

c 

phip(ip)=phiO(ip) 

$  -■0.5*dt/dx*um(ip)*(  phiin(ip2)“phim(ipl)  ) 

$  -0.5*dt/dy*vm(ip)*(  phim( jp2)-phim( jpl)  ) 

c 

c  integrate  U 

c 

up(ip)=uO(ip) 

c 

c  integrate  V 

c 

vp(ip)=vO(ip) 
end  do 

end  do 
c 

c  Apply  the  Periodic  Boundary  Conditions 
c 

do  j=l,hy 

il=node(l , j ) 
i2=node(nx , j ) 
phip(il)=phip(i2) 
up(il)=up(i2) 
vp(il)=vp(i2) 
end  do 

do  i-l,nx 
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il=node(i, 1) 
i2=node(i ,ny) 
phip(il)=phip(i2) 
up(il)=up(i2) 
vp(il)=vp(i2) 
end  do 

1000  continue 

return 

end 

- - * 

♦This  subroutine  writes  the  output.  It  is  currently  set  only  to 
♦print  the  concentration  (or  color)  function  at  each  node  point. 
♦Written  by  F.X.  Giraldo  on  2/95 

- - - 

subrout ine  output (phi , u , v , npo in , t ime , nx , ny ,mxpo i ) 

implicit  real*8(a-h,o-z) 

dimension  phi(mxpoi),  u(mxpoi),  v(mxpoi) 

pi=4. 0*atcai(l  .0) 

writed,  '  (2(i6,lx)  ,el6.8)  ' )  ,nx,ny , time/(2 . 0+pi) 
writed,  '  (el2.4)  d(phi(ip) ,  ip=l,npoin) 
return 
end 


♦These  next  4  subroutines  construct  the  Hermitian  Interpolation  functions 
♦required  by  the  semi-Lagrangian  method. 

♦Obtained  from  Numerical  Recipes. 

♦Written  by  F.X.  Giraldo  on  2/95 

- - 

subrout ine  spl ie2 ( df , f . coord , nx , ny , node , mxpo i , imax ) 
implicit  real*8(a“h,o-z) 
parameter  (nmax=3000) 
c 

c  global  arrays 

c 

dimension  coord(mxpoi ,2) ,  f(mxpoi),  df(mxpoi) 
integer  node ( imax, iniax) 
c 

c  local  arrays 

c 

dimension  x(nmax),  y(nmax) ,  y2(nmax) 
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do  j=l,ny 
do  i=l,nx 

y(i)=f (node(i, j)) 
x(i)=coord(node(i , j ) ,  1) 
end  do 

call  spline(y2 ,y ,x,nx, idiiin,nmax) 
do  i=l,nx 

df (node(i, j))=y2(i) 
end  do 
end  do 

return 

end 

- - * 

subroutine  spl in2 (f d , coord , f , df , xd , yd , nx , ny , node , mxpo i , imax ) 
implicit  real*8(a-h, o-z) 
parameter  (nmax=3000) 
c 

c  global  arrays 

c 

dimension  coord(mxpoi ,2) ,  f(mxpoi),  df(mxpoi) 
integer  node ( imax , imax) 
c 

c  local  arrays 

c 

dimension  x(nmax) ,  y(nmax),  y2(nmax),  y22(nmax) 

do  j=l,ny 
do  i=l,nx 

y(i)=f (node(i, j)) 
y2(i)=df (node(i, j ) ) 
x(i)=coord(node(i, j) ,1) 
end  do 

call  splint (fd,x,y,y2,xd,nx,nmax) 
y22(j)=fd 
end  do  ' 
do  j-l,ny 

x( j )=coord(node(l , j ) ,2) 
y(j)=y22(j) 
end  do 

call  spline(y2,y,x,ny,idum,nmax) 
call  splint (f d , x , y , y2 , yd , nx ,nmax) 

return 
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end 


subrout ine  spline ( y 2 , y , x , n , idum , mx ) 
implicit  real*8(a-h, O'-z) 
parameter  (nmax=3000) 

dimension  y(mx) ,  y2(rax),  x(mx),  u(nmax) 
if  (mx . gt . nmax )  then 

write(*,'("  Must  expand  NMAX  in  Subroutine  Spline")^ 
stop 
endif 

y2(l)=0.0 
u(l)=0.0 
do  i=2,n-l 

sig=(x(i)-x(i-l) )/ (x(i+l)“x(i-l) ) 
p=sig*y2(i-l)  +  2.0 
y2(i)=(sig  -  1.0)/p 

u(i)=(6.0*((y(i+l)“y(i) )/(x(i+l)”x(i) )  ”  (y(i)-y(i-'l) ) 

$  /(x(i)-x(i-l)))/(x(i+l)-x(i~l))-sig*u(i-l))/p 

end  do 
y2(n)=0.0 
do  i=n-l,l,-l 

y2(i)=y2(i)*y2(i+l)  +  u(i) 
end  do 

return 

end 

- - - 

subrout ine  splint (yd , x , y , y2 , xd , n ,mx) 
implicit  real*8(a-h, o-z) 
dimension  x(mx) ,  y(mx) ,  y2(mx) 


il=l 

i2=n 

10  if  (i2-il.gt.l)  then 
im=(i2+il)/2 
if  (xd.gt .x(im))  then 
il=im 
else 

i2=im 
endif 
goto  10 
endif 
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xl=x(il) 

x2=x(i2) 

dx=x2-xl 

a=(x2-xd)/dx 

b=(xd-xl)/dx 

yd=a*y(il)  +  b*y(i2) 

$  +  ((a**3‘-a)*y2(il)  +  (b+*3-b)*y2(i2) )*(dx**2)/6.0 

return 
end 

♦ - ♦ 

♦This  subroutine  updates  the  arrays  PHIM,UM, VM, ALFM,  PHI0,U0,V0, ALFO 
♦Written  by  F.X.  Giraldo  on  2.95 

♦ - * 

subrout ine  update (phim , phiO , phip , um , uO , up , vm , vO , vp , alf m , alf 0 , 

$  npoin,mxpoi) 

implicit  real^8(a-h,o-z) 

dimension  phim(mxpoi),  phiO(mxpoi),  phip(mxpoi) 
dimension  um(mxpoi) ,  uO(mxpoi),  up(mxpoi) 
dimension  vm(mxpoi),  vO(mxpoi),  vp(mxpoi) 
dimension  alfm(mxpoi , 2) ,  alf 0(mxpoi,2) 


c 

c 

c 

c 

c 

c 


c 

c 

c 


Loop  through  all  the  nodes  and  update 


do  ip=l,npoin 

Update  F(x“2+ alpha , t-dt )=F(x“alpha, t ) 


ph  im  (  ip )  =ph  iO  ( ip ) 

um(ip)=uO(ip) 

vm(ip)=vO(ip) 

alf m ( ip , 1 ) =alf 0 ( ip , 1 ) 

alf m ( ip , 2 ) =alf 0 ( ip , 2 ) 


Update  F(x-alpha,t )=F(x,t+dt) 


phiO(ip)=phip(ip) 
uO(ip)=up(ip) 
vO(ip)=vp(ip) 
end  do 
return 
end 
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A  Matlab  Program  to  plot  the  output  showing  the  cone  at  various  time  is 
given 

clg 

i=0; 

j=0; 

f id=f open( ^matlab . out ’ , ’r O ; 

ab=f scanf  (f  id,  ^'/,d\  1) ; 

isets=ab(l) 

while  (i  <  isets) 

ab=f  s  canf  (fid,'  7,dy,d'/.f ' ,  3 ) ; 

nx=ab(l) 

ny=ab(2) 

hour=ab(3) 

count=nx*ny ; 

jj=i: 

while  (j  <  ny) 

[aa,  count]=f  scanf  (fid,  ’  */,e'/,e*/,e'/,e*/,e  ’  ,nx)  ; 
ab(j j : j j+count-l)=aa; 

jj=j j+count; 
end; 

u=reshape(ab,nx,ny) ; 
v=u '  ; 

c=contour(v, 10) ; 
clabel(c); 

title([' concentration  after  ' ,num2str(hour) , '  revolutions']) 
print  c  -dps  -append 
i=i+l ; 
j=0; 
end; 
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Conclusion 

We  have  developed  a  bilinear  finite  element  Fortran  code  to  solve  the  two 
dimensional  advection  equation  on  a  unix-based  SUN  Sparc  10  workstation.  The 
stability  of  the  method  is  analyzed.  We  have  also  developed  a  semi-Lagrangian 
finite  element  code.  These  codes  were  experimented  with  in  solving  the  cone 
test  problem.  It  is  clear  from  the  plots  that  the  semi-Lagrangian  is  superior  to 
the  Eulerian  finite  elements,  since  the  cone  rotates  back  to  its  position  without 
leaving  a  noisy  trail  behind. 
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Appendix 


2 


4 


6- 


concentration  after  0  revolutions 


Explicit 

Tl^  —  —  2 1 

At  =  .0125 
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50 


J _ I - 1 - 1 - 1 - 

4  6  8  10  12  1 


Explicit 

rix  =  ny  —  21 

Af  =  .0125 


51 


10  12  14  16  18 


Crank-Nicholson 

=  Tiy  =  21 


At 


.0125 


20 
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concentration  after  0.5  revolutions 

- 1 - 1 - 1 - 1 - r 


_^2.84 


J _ I _ I _ L_ 

8  10  12  14 


Crank-Nicholson 

rix  =  fly  =  21 


concentration  after  0  revolutions 


Fully  implicit 

Tl^  -  Uy  —  2  1 

At  =  .0125 
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concentration  afl 


Fully  impli 
n*  —  — 

At  =  .012 
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10  12  14  16  18  20 


Senii-Lagrangian 

tlx  =  Uy  =21 

At  =  .0125 
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8  10  12  14  16  1 


Semi-Lagrangian 

Tl^  *“  — —  2  ] 

At  =  .0125 


8  20 
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10 


12 


Sem  i-  Lagrangi  an 

nj,  =  rij,  =  2 1 


At  =  .0125 
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